Many-body instability of Coulomb interacting bilayer graphene: RG approach 
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Low-energy electronic structure ol (unbiased) bilayer graphene is made ol two Fermi points with 
quadratic dispersions, if trigonal-warping and other high order contributions are ignored. We show 
that as a result of this qualitative difference from single-layer graphene, short-range (or screened 
Coulomb) interactions are marginally relevant. We use renormalization group to study their effects 
on low-energy properties of the system, and show that the two quadratic Fermi points spontaneously 
split into four Dirac points, at zero temperature. This results in a nematic state that spontaneously 
breaks the six-fold lattice rotation symmetry (combined with layer permutation) down to a two-fold 
one, with a finite transition temperature. Critical properties of the transition and effects of trigonal 
warping are also discussed. 



The ability to predict the nature of the low tempera- 
ture state of an interacting quantum system is one of the 
main goals of condensed matter theory. Nevertheless, 
despite ongoing effort, no single method has proved uni- 
versally sufficient and experimental input is essentially 
inevitable. 

Under special circumstances, however, progress can be 
made. In particular, in non-interacting systems with 
susceptibilities diverging as the temperature approaches 
zero, the inclusion of arbitrarily small interaction can be 
shown to lead to a finite, but also arbitrarily small tran- 
sition temperature. The method of choice in this case is 
the renormalization group (RG), which has the virtue of 
unbiased determination of the leading instability^. 

In this paper we apply the RG method to the bilayer 
graphene with Bernal stacking^''^'**^ . While in general, 
the motion of the non-interacting electrons in such po- 
tential does not lead to diverging susceptibilities since the 
energy spectrum has two sets of four Dirac points in the 
corners of the Brillouin zone (due to trigonal warping 
if only nearest neighbor hopping is considered each set of 
four Dirac points merges into a single degenerate point 
with parabolic dispersion (See Fig. [1]). As the nearest 
neighbor hopping amplitudes are the largest, the latter 
is the natural starting point of theoretical analysis^ii. 
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We start with the tight-binding Hamiltonian for elec- 
trons hopping on the bilayer honeycomb lattice with 
Bernal stacking 



FIG. 1: (Upper left) Honeycomb bilayer unit cell. Atoms 
in the lower layer (2) are marked as empty (black) circles, 
atoms in the upper layer (1) are filled (red) circles. As a 
starting point, only the intralayer nearest neighbor hopping 
amplitudes t and the interlayer hopping amplitudes t± are 
considered. (Upper right) Constant energy contours of the re- 
sulting dispersion, with minima aX K = y and K' points 
and maximum at P point. (Lower left) The energy dispersion 
of the four bands along the vertical cut in the Brillouin zone. 
The band splitting at the K (and K') points is t±. (Lower 
right) Magnification of the dispersion (in units of t) near the 
degeneracy point (solid black) as well as the dispersion in the 
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where, in the nearest neighbor approximation, the (real) 
hopping amplitudes t connect the in-plane nearest neigh- 
bor sites belonging to different sublattices and, for one 
of the sublattices, also the sites vertically above it with 
amplitude t±. Since there are four sites in the unit cell, 
there are four bands whose dispersion for the above model 



comes from the solution of the eigenvalue problem: 
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FIG. 2: Diagrams appearing at 1-loop RG. The vertices are 
either 6^/3 or E^g. 



We find E{k) = ± ± ^MkP + iti) 

2 cos ^ 



with dk = 
Two of the bands are 



gapped (at K, K' by ij^) and become separated from the 
low energy pair which touches at k = (See FiglT]). The 
resulting density of states at zero energy is therefore fi- 
nite. 

The repulsive interaction V{r — r') in Eq.([T]) is taken to 
have a finite range f which is however much larger than 
the lattice spacing a. This is assumed to be the cor- 
rect starting point, since the full Coulomb interactions is 
screened^, at low energy due to the finite density of states. 
The analysis starting from the l/|r — r'| interaction will 
be postponed to a future publication. 

Following Nilsson et. al^ we project out the gapped 
bands. The resulting low energy effective (imaginary 
time) action (which includes both K and K' valleys) is 



S = 



drd^r 
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where the four component Fermi (Grassman) fields 
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The Pauli matrices CTj act on the layer indices 1-2 and the 
T matrices act on the valley indices K-K'. The effective 
mass is m = 2ij_/(9t^), and represents -y— copies of 
the four component pseudo-spinor. = 4 for spin 1/2, 
and e.g. for s = 1, . . . iV, ip^T,''ip{r,T) = ^pl^T,^p'ipi3s. 
Note that S's have the same multiplication table as the 




FIG. 3: RG flow diagram of the ratios gi/ga and 52/(73 for 
(73 < 0. While the ratio gi/ga flows to zero (even if the start- 
ing point is (?2 = 33 = and gi / 0), the ratio 52/(73 flows 
to a fixed value, indicating two stable and one unstable rays 
with slopes mi « —0.525, 7713 ~ 13.98 and m2 « 0.545, re- 
spectively. 



Pauh cr's: S'^S'' = US^ 



ic^^aS and are traceless, 
too. A is a momentum cutoff which restricts the modes 
to the vicinity of the K-K' points and whose order of 
magnitude is < ^/2'mt±. 

The coupling constant gi = J d'^rV{r), i.e. it is the 
q = Fourier component of V{r). The coupling con- 
stants 52 and gs are zero in the starting action, but as 
will be shown next, they get generated in the momentum- 
shell RG^, and therefore they are made explicit in the 
original action. 

From simple power-counting, the (engineering) scal- 
ing dimension of the field ip is L"^ and for r. This 
makes gi, 52 and 53 marginal (at the tree- level) and the 
question is how they flow upon inclusion of the loop cor- 
rections. To answer this we note that all possible Wick 
contractionsi of four-fermion operators correspond to the 
diagrams in the Figure 0. The RG equations obtained 
by integrating fermion modes within a thin shell A and 
A/s (centered at the K point), and are: 
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[-A{N - l)gl + 4.g| + 4.gi.g2 - I25253] ^(8) 



= [-(51-53)' 
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While the above equations cannot be solved in a closed 
form, it is possible to fully analyze the qualitative nature 
of the RG flows. Such analysis is facilitated by the ob- 
servation that 
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TABLE I: The susceptibility coefficients A,B,C in Eg. dTTll 
for different particle- hole order parameters ijj''Oiip. In the 
physical case A'^ = 4. 



TABLE IL The susceptibility coefficients A', B' , C in Eq.(fT8| 
for different particle-particle order parameters iJaa-O^^iipfja-' ■ 



which means that, unless gi — g2 — 93 = when the 
equality holds, 53 strictly decreases under RG rescaling. 
We can therefore trade the parametric dependence on s 
of gi and 52 for their dependence on 53 and retain the 
direction of the RG flow. For 53 < (> 0), an increase in 
d log s therefore corresponds to an increase (decrease) in 
Since the system is autonomous, we can eliminate 
log s and arrive at a system 

^ ^ f (E1 92 

dgs ygs' ga 
dm ^ fgigi 
dg3 ygs' gs 
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The system of Eas. (fTO|) -([TT |) is in turn homogeneous 
and can therefore be written as 



53 



93 



d^ 
dgs 

93 

dg3 



gi 
93 



91 52 
gs ' gs 



.92 , gi 92 
— + .9 — , — 
.93 V53 53 



(14) 
(15) 



The above system has three fixed points, all of which 
have 51/53 = 0, while 92/93 = mi,TO2,m3. As shown 
in the Fig. ([3]), mi « —0.525 and 1713 « 13.98 are sinks, 
while 7712 ~ 0.545 has one attractive direction and one 
repulsive. This means that once 93 gets to be negative, 
only g2 and 93 become important (their ratio being fixed) 
while gi is too small compared to 93. To see that this is 
indeed what happens if the starting point is 171(3 = 1) > 
and 92(3 = 1) = 93(3 = 1) = 0, note that the Eqs.(I7]- 
[9]) imply that finite gi generates finite and negative 93 
upon first iteration while 52 remains zero until the second 
iteration. This means that we start with gi/ga — > — cxd 
and 32/53 = which is below the (red) separatrix, thus 
the flow is into the region of attraction of mi (Fig. ([3])). 

From Eqs.(17][ni) we see for the fixed ratios 51/53 — 
and 52/53 = ruj, 53 becomes large and negative, indi- 
cating a runaway flow. Given the flow of the coupling 
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FIG. 4; Numerical integration of the susceptibilities in 
Eq.([T7ll for gi{s = 1) = 0.01 and gzis = 1) = gais = 1) = 0. 
The strongest divergence is towards the nematic order. (In- 
set) Numerically determined nematic transition temperature 
in units of cutoff Ta < as a function of the dimensionless 
coupling gi^. 



constants we can determine the susceptibilities towards 
the formation of ordered states. In particular, we con- 
sider coupling the fermions to external sources, which 
correspond to the possible broken symmetry states. We 
therefore have additional terms in the action: 

AS = -A^^ J dTd\ij^O,i^{r,T) 

- A^p' / drd^riP^^Ol^f^i^p^'ir^T) (16) 



Such terms, with infinitesimal A's explicitly break the 
symmetry and so are relevant operators. The question 
of instability is answered by finding the renormalization 
of the vertices^°. The one which diverges first deter- 
mines the broken symmetry states. After a straigthfor- 
ward calculation we find that for a general particle-hole 



order parameter Oi 

TQ = (To = 1, 
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where the coefficients A, B, and C are given in the Table 
m Similarly, for a general particle-particle order parame- 



ter %pc 



pp^ren '—^pp 



l + [A'gi + B'g2 + C'g3]—liis 



(18) 



where the coefficients A', B', and C are given in the 
Table m 

The instabihty towards a particular order occurs at an 
energy scale (i.e. temperature) at which the correspond- 
ing coefficient of the Ins in Eas. (|17m"8|) diverges. Since 
TV = 4 and the fixed point value of ^2/33 ~ —0.525, with 
gs large and negative, it can be seen from Table H] that 
the instability appears in the "E^^v channel, which as we 
discuss next corresponds to a nematic order. The numer- 
ical integration of the RG equations ((TE]) starting with 
gi{s = 1) > and g2{s = 1) = (73(5 = 1) = shown in 
Fig.(|4]) indeed confirms that the susceptibility diverges 
fastest in this channel. Within the continuum model and 
in weak coupling, the instability is therefore towards the 
order parameter, which we can parametrize by a complex 
field 

A„e™(r) = A,(r) + tAy{r) = {^Hr) (S^ + iS^) ^(r)). 

To see that this is indeed a nematic order, note that at 
q = (1) it is translationally invariant and (2) even under 
rotations by tt. In fact, as the low energy Hamitonian is 
invariant under arbitrary rotations by an angle a, i.e. 
W{a)nU{a) = n, where = e''"^- e''"'^', = 



x-§- — v-i-, we find that under a rotation by a 

ay ^ ax ' ^ 



i(r) A„e,„(r)e^ 



This shows that the order parameter is even under rota- 
tions by TT and odd under rotations by 7r/2, which makes 
it nematic. For uniform A„em(r) the quadratic degener- 
acy point is split into two (massless) Dirac points by an 
amount proportional to the magnitude of the order pa- 
rameter and the direction given by the nematic director. 

The presence of the underlaying lattice further breaks 
the full rotational symmetry of the long distance effective 
Hamiltonian down to hexagonal symmetry centered on 
0-2 — ^1 site, where the standard operations of Cgu must 
be accompanied by the appropriate layer permutations. 
The two components of the order parameter, which give 
finite expectation values of, for instance, Aj;(r) = 



and 




s=± J 



+ h.c. 



form a two dimensional representation of the hexagonal 
group. Note that the nematic order parameter remains 
even under 7r-rotation followed by the layer permutation. 

From the arguments above we expect that the lat- 
tice has an important effect on the critical nature of the 
phase transition, which would otherwise be of Kosterlitz- 
Thousless kind. The reason is the existence of the third 
order invariant Ai^ — 3Aa;A^. As a result the finite tem- 
perature phase transition should be described by the ef- 
fective Hamiltonian 



Unem = -'-'^ COS [2(6* (x) - 6l(x'))] + /i ^ cos[66'(x)];21) 



(xx'> 



where A^(x) -h jAy(x) = e^*^^^), 9 £ (0,27r] and the 
sum runs over the vertices of the triangular sub-lattice 
spanned by ai sites. This corresponds to the p = i case of 
the two dimensional planar model studied by Jose et.al.— 
and the concomitant absence of the Gaussian spin-wave 
phase. Instead there is a continuous transition between 
the low temperature phase where the director locks into 
one of three values and a high temperature phase where 
vortices unbind. Such transition is believed to belong to 
the 2D three-state Potts model universality classic with 



exponent: 



,13 



5/6 and 77 = 4/15. 



Finally, we discuss the effects of the trigonal warping 
which splits each of the quadratic degeneracies into four 
massless Dirac points, which were ignored up to now. If 
we denote the energy scale associated with such terms as 
Ttrig, below which the dispersion must be modified, then 
the transition will still occur provided that the mean- 
field transition temperature Tc estimated from the above 
model and plotted in the inset of Fig. ([3]) satisfies ^ 
Ttrig. For screened Coulomb interactions- 5137 ~ C(l): 
leading to < tj^. Since the current estimates of Ttrig 
are of the same order of magnitude'^'*, the ultimate test 
is experimental. 
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